SoSe2021
iris lm()Wir modelliere die Kronblattlänge als Funktion der Kelchblattlänge:
mod <- lm(formula = Petal.Length ~ Sepal.Length, data = iris) mod
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Coefficients:
(Intercept) Sepal.Length
-7.101 1.858
coef(mod)
(Intercept) Sepal.Length -7.101443 1.858433
Wie lautet die vollständige lineare Regressionsgleichung?
\[y = a + b*x\]
Die Stichprobenkennwerte \(a\) und \(b\) repräsentieren die Punktschätzer der Regressionsparameter \(\alpha\) und \(\beta\)!
iriscoef(mod)
(Intercept) Sepal.Length -7.101443 1.858433
Für die beobachteten Stichprobenwerte gilt:
\[Petal.Length_i = -7.1 + 1.85*Sepal.Length_i + \epsilon_{i},~~\text{wobei}~~\epsilon_i ~N(\mu, \sigma^2)\]
summary() Funktionsummary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Wo ist a? Wo ist b? Und was bedeuten die anderen Zahlen…?
Stelle Dir vor, jeder Datenpunkt repräsentiert einen individuellen Fisch und alle Datenpunkte zusammen repräsentieren eine Population.
Wenn Du jetzt mehrere Zufallsproben von 10 Fischen nimmst, wie ähnlich wären sich dann die geschätzten Koeffizienten?
Rate: Welche sind sich viel ähnlicher und welche sind näher an den wahren Populationskoeffizienten dran?
‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41
Lass uns jetzt die Probengröße auf 100 Fische steigern!
‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41
Wie ähnlich sind sich die geschätzten Koeffizienten bei schwacher Streuung und einem Stichprobenumfang von 10 Fischen?
‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5
Lass uns jetzt die Probengröße auf 100 Fische steigern!
‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5
Die Unbestimmtheit der geschätzten Steigung und Achsenabschnitt
\[s_\bar{a}=\sqrt{\frac{MS_{Residual}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}\sum x_i^2}{n*\sum(x_i-\bar{x})^2)}}\]
\[s_\bar{b}=\sqrt{\frac{MS_{Residual}}{SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}}{SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}}{\sum(x_i-\bar{x})^2}}\]
Die Gesamtstreuung (\(SS_Y\) oder \(SS_{Gesamt}\)) kann in 2 Komponenten zerlegt werden:
erklärbare Streuung = \(SS_{Regression}\) + nicht erklärbare Streuung = \(SS_{Residuen}\)
summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16→ Punktschätzung für den Achsenabschnitt: -7.10144 ± 0.50666
→ Punktschätzung für die Steigung: 1.85843 ± 0.08586$
confint()a <- -7.10144
b <- 1.85843
ci_a <- 0.50666 * qt(p = 0.975, df=148)
ci_b <- 0.08586 * qt(p = 0.975, df=148)
mat <- matrix(
c(a-ci_a, a+ci_a, b-ci_b, b+ci_b),
nrow = 2, byrow = TRUE)
rownames(mat) <- c("(Intercept)", "Sepal.Length")
colnames(mat) <- c("2.5%", "97.5%")
mat
2.5% 97.5% (Intercept) -8.102662 -6.100218 Sepal.Length 1.688760 2.028100
confint(mod, level = 0.95)
2.5 % 97.5 % (Intercept) -8.102670 -6.100217 Sepal.Length 1.688772 2.028094
(t_a <- -7.10144 / 0.50666)
[1] -14.01618
(t_b <- 1.85843 / 0.08586)
[1] 21.64489
pt(t_a, df = 148)
[1] 3.065741e-29
pt(t_b, df = 148, lower.tail = FALSE)
[1] 5.224047e-48
summary() outputsummary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Wenn ein p-Wert sehr niedrig ist ( < 0.05), bedeutet das, dass
Die statistische Kenngröße \(R^2\) ist ein Maß für den Anteil der Variabilität in y, welcher durch das lineare Modell erklärt werden kann:
\[R^2 = \frac{SS_{Regression}}{SS_{Gesamt}} = \frac{\sum(\hat{y}_i-\bar{y})^2}{\sum(y_i-\bar{y})^2}\] \[1 \geq R^2 \geq 0\]
summary() outputsummary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)
Residuals:
Min 1Q Median 3Q Max
-2.47747 -0.59072 -0.00668 0.60484 2.49512
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.10144 0.50666 -14.02 <2e-16 ***
Sepal.Length 1.85843 0.08586 21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared: 0.76, Adjusted R-squared: 0.7583
F-statistic: 468.6 on 1 and 148 DF, p-value: < 2.2e-16Es wird angenommen, dass
p
Angenommen Du hast 12 Messwerte für die Variable x und y und willst den (positiven) Einfluss von x quantifizieren:
Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:
Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:
Führe eine Modellvalidierung durch und erstelle die sog. Diagnostikplots mit den Residuen:
→ Es ist einfacher den Ausschluss beeinflussender Datenunkte zu rechtfertigen, wenn eine große Cook’s Distanz und eine große Leverage für den Datenpunkt vorliegt.
p
par(mfrow = c(2,2)) plot(mod)
Du kannst beide Arten von Residuen berechnen mit
residuals(model) (oder resid()) vom ‘stats’ Paket → ergibt einen Vektor mit den gewöhnlichen Residuen
rstandard(model) vom ‘stats’ Paket → ergibt einen Vektor mit den standardisierten Residuen
add_residuals(data, model, var = "resid") vom ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘resid’ mit den gewöhnlichen Residuen zu deinem ‘data frame’ hinzu; hilfreich beim ‘piping’ von Funktionsaufrufen!
Du kannst Dir die vorhergesagten Werte für dein Modell über 3 Funktionen ausgeben lassen:
fit(model) → gibt einen Vektor mit den angepassten Werte zurück, die vom Modell für die Werte der erklärenden Variablen vorhergesagt werdenpredict(model) von dem internen ‘stats’ Paket → verwendet Informationen aus dem angepassten Modell, um eine Glättungsfunktion für die Darstellung einer Linie im Streudiagramm zu erstellen. Der Funktion kann ein data frame mit Sequenzen von Werten aller X Vaiablen, die im Modell sind, übergeben werden (Argument newdata).add_predictions(data, model, var = "pred") von dem ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘pred’ mit den vorhergesagten Werten zu deinem Datensatz hinzu; hilfreich bei ‘piping’ Operationen!iris %>% modelr::add_predictions(mod) %>% ggplot(aes(x = pred, y = Petal.Length)) + geom_point() + geom_abline()
iris2 <- iris %>% modelr::add_residuals(mod) p <- ggplot() + geom_hline(yintercept = 0, color = "red") p1 <- p + geom_point(aes(x = Sepal.Length, y = resid), iris2) p2 <- p + geom_point(aes(x = Sepal.Width, y = resid), iris2) p3 <- p + geom_point(aes(x = Petal.Width, y = resid), iris2) gridExtra::grid.arrange(p1,p2,p3, nrow = 1)
→ Die Residuen sollten zufälligt verstreut sein. Ein Muster deutet auf eine Fehlspezifikation des Modells hin, wie hier zum Beispiel.
Verwende geom_smooth und setze das Methodenargument als Name der Modellierungsfunktion (hier: ‘lm’). Falls Du den Standardfehler nicht anzeigen willst, setze zusätzlich se = FALSE:
iris %>% ggplot(aes(x = Sepal.Length, y = Petal.Length)) + geom_point() + geom_smooth(method = "lm", se = FALSE)
iris %>% modelr::add_predictions(mod) %>% ggplot(aes(x = Sepal.Length)) + geom_point(aes(y = Petal.Length)) + geom_line(aes(y = pred))
Erkunde die Shiny App auf der folgenden Folie und beantworte diese zwei Fragen:
Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)